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Abstract 



The ground state geometries of some small clusters have been obtained 
via ab initio molecular dynamical simulations by employing density based en- 
ergy functionals. The approximate kinetic energy functionals that have been 
employed are the standard Thomas-Fermi {Tj<f) along with the Weizsacker 
correction T\y and a combination F(N £ )Ttf + Tyy. It is shown that the 
functional involving F(N e ) gives superior charge densities and bondlengths 
over the standard functional. Apart from dimers and trimers of Na, Mg, Al, 
Li, Si, equilibrium geometries for Li n Al,n = 1,8 and Al\^ clusters have also 
been reported. For all the clusters investigated, the method yields the ground 
state geometries with the correct symmetries with bondlengths within 5% 
when compared with the corresponding results obtained via full orbital based 
Kohn-Sham method. The method is fast and a promising one to study the 
ground state geometries of large clusters. 
PACS Numbers : 71.10, 31.20G, 02.70N, 36.40 
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I. INTRODUCTION 



During the last few years the technique of first principles molecular dynamics (MD), 
initiated by Car and Parrinello (CP) has emerged as a powerful tool for investigations 
of structural, electronic and thermodynamic properties of large scale systems. The standard 
implementation of this method which is based on density functional theory is via Kohn-Sham 
(KS) orbitals. Such orbital based algorithms scale as N%, N a being the number of atoms in 
the system. Quite clearly, such methods turn out to be computationally expensive for system 
sizes over about 100 atoms M. Recently, approaches based on total energy functionals, which 
depend on charge density only or orbital free density functionals have been proposed . 
These methods are based on approximate representation of kinetic energy (KE) functionals 
and offer an attractive alternative for investigating large scale systems. Since the method is 
orbital free i.e there are no wavefunctions to handle, there is no computationally expensive 
orthogonality constraint and the methods scale linearly with system size. In addition, these 
methods are shown to yield stable dynamics even with large timesteps, a highly desirable 
feature for molecular dynamics simulations. 

It is clear that the utility of these methods is critically dependent on their ability to 
investigate the systems of interest with acceptable accuracy, at least for a class of physical 
properties. Madden and coworkers have investigated structural and thermodynamic prop- 
erties of some simple metals with considerable success. For example, the dynamic structure 
factor of liquid Sodium and static structure factor, vacancy formation energy, free energies 
of point defects as well as phonon dispersion curves of Sodium j|,[7| are well described by 
this method. The method has also been applied for ground state configurations of c-Si and 
H/Si (1 0) surface || and for geometries of some silicon clusters |§ and a good agreement 
has been found with experiments as well as with other calculations. However majority of 
the calculations reported so far have been performed on extended systems. 

In the present work, we focus our attention on studying the ground state and ener- 
getically low lying structures of clusters, a field of current interest. Obviously, due to the 
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approximate nature of the KE functionals the bondlengths and binding energies will not 
be obtained with the same level of accuracy as the KS orbital based methods. However, 
it is of considerable interest to examine whether such a method is capable of yielding the 
correct shapes (i.e the right symmetries) of clusters by employing Car-Parrinello simulated 
annealing methods. If desired the KS method can then be used to search the local minimum 
around structures obtained by Orbital Free Method (OFM) in 'quenching' mode. This can 
be a computationally tractable way to avoid the long and costly simulated annealing runs 
of the orbital based KS molecular dynamics. 

Towards this end we have carried out a number of calculations on a variety of represen- 
tative small clusters of simple metals. Specifically, we have investigated dimers and trimers 
of Na, Mg, Al, Li, Si, small clusters of Na n , (n = 6,8), Li n Al, (n = 1,8) and Ali 3 . These 
systems are representative of the small metal atom clusters of current interest and more 
accurate KS based results have been reported. Hence, it is possible to make an assessment 
of the present method by comparing the bondlengths and geometries with the reported ones. 

The question of appropriate choice of kinetic energy functionals has been addressed by 
Smargiassi and Madden || . They have investigated a family of kinetic energy functionals 
which incorporate exact linear response properties. All such KE functionals are based on 
the Thomas Fermi (TF) and the Weizsacker correction term. Since our interest is in finite 
size systems, we have chosen to use simple KE functionals. These functionals have been 
previously used in the study of atoms and molecules. However, it must be mentioned 
that significant progress has been made towards improving the KE functionals notably by 



DePristo and Kress, Wang and Teter 10.11 



In the next section we briefly discuss the method, the KE functionals used and give the 
relevant numerical details. This is then followed by the results and discussion. 



II. FORMALISM AND COMPUTATIONAL DETAILS 
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A. Total Energy Calculation 



The total energy of a system of N e interacting electrons and N a atoms, according to the 
Hohenberg-Kohn theorem [12|,13|], can be uniquely expressed as a functional of the electron 
density p(r) under an external potential due to the nuclear charges at coordinates R n , 



E 



p, {R n }] = T[p] + E xc [p] + E c [p] + E ext [p, {R n }] + ^({R n }), 



(1) 



where E xc is the exchange-correlation energy, E c is the electron-electron Coulomb interaction 
energy. The electron-ion interaction energy E ext is given by 



E 



est 



p,{R n } = / V(r)p(r)d 3 r 



(2) 



where ^(r) is the external potential, usually taken to be a convenient pseudopotential |T4| 



The last term in Eq. (1), En, denotes the ion-ion interaction energy. The first term in Eq. 
(1), the KE functional, is usually approximated as 



T[p] = T TF [p}+T w [p\ 



(3) 



where Ttf[p\ is the Thomas- Fermi term, exact in the limit of homogeneous density, and has 
the form 



T TF [p] = 1(37T 2 )I / p(rfW 



(4) 



and T w [p] is the gradient correction due to Weizsacker, given as 

A r Vp(r) • Vp(r)d 3 r 



Twlfi] 



p(r) 



(5) 



which is believed to be the correct asymptotic behavior of T[p] for rapidly varying densities. 
Instead of A = 1, the original Weizsacker value, A = | and A = g are also commonly 
used. It has been argued that for rapidly varying densities, which is the case for finite size 
clusters a more appropriate kinetic energy functional would be a following combination of 



these two terms 16 
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T[p] = F(N e )T TF [p}+T w [p\ (6) 

where the factor F(N e ) [17] is 

f(au = (i-^)(i-4+4) m 

V 7V e /V jus/ 



with optimized parameter values Ai = 1.314 and A<i = 0.0021 ||18|| . This functional which 
includes the full contribution of the Weizsacker correction describes the response properties 
of the electron gas well. This functional has been used for investigating atoms and molecules 
with reasonable success. 

We briefly describe our procedure, details of which can be found in ||. The total energy 
of the system (Eq. (1)), is minimized for fixed ionic positions using the conjugate gradient 



method |19| which forms the starting point for molecular dynamics. The trajectories of 
ions and the fictitious electron dynamics are then simulated using Lagrange's equations of 
motion which are solved by Verlet algorithm JT|. The stability of CP dynamics has been 
discussed in Q in the context of density based methods and timesteps of the order to 50 
a.u. have been successfully used. We have verified that by appropriate adjustment of the 
fictitious electron mass the CP dynamics remains very stable for over 10000 iterations with 
a timestep of 40 a.u. in the present calculations of clusters. Typically, for free dynamics the 
grand total energy which is the sum of the kinetic energy of ions, kinetic energy of electrons 
and the potential energy of the system remains constant to within 10~ 5 a.u. 

For the calculations of the ground state structures for dimers and trimers of Na, Mg, 
Al, Li and Si a periodically repeated unit cell of length 26 a.u. with a 54 x 54 x 54 
mesh and timestep At ~ 10 to 20 a.u. was used. For the rest of the small clusters the 
calculations were done on a unit cell of length of 30 a.u with a 54 x 54 x 54 mesh. We 
have chosen to use the plane wave expansion on the entire fast fourier transform mesh 
without any truncation yielding the energy cutoff of 95 Rydberg. It must be mentioned 
that, due to the orbital free calculations the number of fast fourier transforms per iteration 
are constant irrespective of the number of electrons in the system. For clusters, the ground 
state configurations are obtained either by starting with different initial configurations and 



then quenching the structures or by dynamical simulated annealing where the cluster is 
heated to 300 — 350° K and then cooled very slowly. In all the cases the stability of the final 
ground state configurations has been tested by reheating the clusters and allowing them to 
span the configuration space for a few thousand iterations and then cooling them to get the 
low energy configurations. 

III. RESULTS AND DISCUSSION 

In this section, we first discuss the results for the equilibrium bondlengths and binding 
energies of dimers and trimers of Na, Mg, Al, Li, Si along with their KS results. All the 
results presented here are obtained with energy convergence up to 10~ 13 for total energy 
minimization. 

Table I shows the equilibrium bondlengths and binding energies for dimer and trimer 
systems using different kinetic energy functionals. These results have been compared with 
full nonlocal pseudopotential KS calculations. A few representative results using A = | have 
been given. It can be seen that for A = | the trend is similar to the A = | functional 
and there is no significant improvement in the results. Clearly, the results involving F(N e ) 
functional show significant improvement over A = | (with the exception of Mg) and are 
in reasonable agreement with the bondlengths obtained by the KS method. The error in 
the bondlengths is around 10%. It is known that such methods based on approximate KE 
functionals are not expected to give accurate binding energies. One notable feature of the 
binding energy comparison is the considerable improvement by F(N e ) over A = | (excepting 
again the case of Mg). The results for the Na, Li, Si trimer binding energies are not given 
because these are Jahn- Teller distorted isoscales triangles and the present method yields 
equilateral triangle geometries. Clearly, such density based methods are unable to reach the 
Jahn- Teller distorted geometries. 

The quality of the charge densities obtained by this method can be gauged by comparing 
them with the KS charge density. In Fig. 1 we have plotted the self-consistent charge 
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densities obtained using the functional involving F(N e ) (curve a) and A = | (curve b) with 
the KS charge density (curve c) for Al dimer along the axis joining the atoms. The ionic 
positions are marked by arrows on the plot. The KS charge density has been obtained using 
the identical pseudopotentials and the same cell size as in the case of OFM. Three prominent 
features can be observed. 

1. Overall the F(N e ) functional densities compare very well with the KS densities ex- 
cept at the origin where both the F(N e ) and A = | self-consistent densities show 
overestimation. 

2. At the atomic sites the F(N e ) and KS based densities are very close and nonzero, 
whereas the A = | shows a disturbing feature of almost zero density. 

3. At the peaks on either side of the origin, the KS and F(N e ) charge densities again are 
close, but the charge density by A = | shows considerable overestimation. 

In Fig. 2 we have plotted the superposed free atom charge density (0 th iteration density) 
represented by the curve b and the self-consistent charge density for the functional involving 
A = | by the curve c and F(N e ) by the curve a. The self consistent charge density obtained 
using A = | shows improvement only at the origin. Contrary to this the self consistent 
charge density using F(N e ) shows a significant overall improvement, both at the origin and 
at the peaks on either side of the peak at the origin. To get an idea of the nature of the 
forces obtained by the OFM and KS dynamics, we have given the results for the vibrational 
frequencies for Na, Mg and Li dimers in Table 2. It is gratifying to note that the vibrational 
frequencies obtained by OFM method are in very good agreement with the KS ones. 

To assess the utility and performance of this method, it has been applied to calculate 
the ground state geometries of a range of small clusters. We report here our calculations 
on hetero nuclear clusters of Li n Al,n = 1,8 and a highly symmetric homonuclear cluster of 
Al 13 and clusters of Na n ,n = 6,8 using the F(N e ) functional. The results are compared 
with the ones reported by KS method. 
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The geometries of the heteronuclear Li n Al clusters are shown in Fig. 3 and the 
bondlengths and symmetries in Table III. along with the KS results. Evidently, the present 
method not only reproduces the correct ground state geometry with bondlengths within 5% 



but also reproduces the two key features observed in the more accurate KS calculation [20 



1. The Li n Al clusters for n < 3 are two-dimensional whereas from n > 3 the clusters 
become three-dimensional. 

2. The Al atom gets trapped inside the Li atoms at n = 6. 

It can also be noted that as the number of atoms in the cluster increases, the accuracy 
in the bondlengths appears to improve. However, for the case of Li 3 Al and Li 8 Al we get 
the ground structure configurations to be ideally symmetric rather than slightly Jahn Teller 
distorted geometries of the KS calculation. 

We have also investigated the AI13 cluster since it shows an interesting icosahedral geom- 
etry. The calculations were performed in two different ways. First we started with a highly 
distorted icosahedron and applied the dynamical quenching to get the equilibrium geometry. 
In the second one, we started by placing the Al atoms at the fee lattice points and heated 
the cluster to 300°K and let the system span the configuration space for a few thousand 
iterations. This was then followed by a slow cooling schedule. It is very gratifying to note 
that in both the calculations the correct icosahedron is obtained with a bondlength of 4.88 
a.u. as compared to the KS bondlength of 5.03 a.u. The error in the bondlengths being 
3%. This strengthens our confidence in the ability of the method to reproduce the correct 
ground state geometries with acceptable bondlengths. In addition, we have also obtained 
the ground state geometries for Na6, Nas and Na2o and have verified that the geometries 
obtained are identical to those reported in pif , with the bondlengths differing by about 5%. 
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IV. CONCLUSION 



In this work, we have presented the results obtained by using density based ab initio 
MD for a variety of small clusters and demonstrated that the method using approximate 
KE functionals is capable of yielding bondlengths within an accuracy of 5%. Our calcula- 
tions indicate that the ground state geometries and symmetries of both homonuclear and 
heteronuclear clusters can be obtained within a reasonable accuracy and timesteps of the 
order of 40 a.u. can be used successfully for stable dynamics. The F(N e ) functional is shown 
to give considerable improvement over standard A = | functional both in terms of charge 
densities and bondlengths and is thus recommended. 

We believe the method to be a promising tool in the study of finite temperature and 
dynamical properties of clusters. So far, all the reported OFM calculations have been 
performed using local pseudopotentials only and it would be interesting to implement the 
nonlocal pseudopotentials and study the effect of nonlocality on the bonding and binding 
properties of such clusters. More work is required in this direction and the implementation 
of nonlocality is under consideration. It may be possible to expand the applications of 
OFM by incorporating the nonlocal pseudopotentials and by employing more accurate KE 
functionals. It is hoped that the problems of current interest in the field of clusters like 
fragmentation, dissociation, interaction between clusters, which may involve large number 
of atoms as well as more than one atomic species will be amenable by the present technique. 
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TABLE I. Comparison of the equilibrium bondlengths (in a.u) and binding energies (in 
eV/atom) using the different kinetic energy functionals with the KS self consistent method. 



System Bondlengths Binding energies 

A = § A=§ F(N e ) KS A = § A = ± F(N e ) KS 



Na 2 


5.67 




5.69 


5.66 a 


-0.116 




-0.867 


Na 3 


5.81 


5.99 


5.75 


6.00 6 


-0.207 


-0.281 


-1.286 


Mg 2 


5.79 




4.71 


6.33 c 


-0.195 




-1.432 


Mg 3 


5.94 


5.81 


4.87 


5.93 c 


-0.355 


-0.526 


-2.096 


Al 2 


5.74 




4.14 


4.66 d 


-0.261 




-1.389 


Ah 


5.88 


5.57 


4.32 


4.74 d 


-0.483 


-0.733 


-2.074 


Li 2 


5.87 




5.51 


5.15 e 


-0.102 




-0.891 


Li 3 


6.03 


6.11 


5.58 


5.3 b 


-0.182 


-0.256 


-1.311 


Si 2 


5.35 




3.74 


4.29-f 


-0.371 




-0.56 


Si 3 


5.50 




3.92 


4.10 / 


-0.651 




-0.938 



-0.71" 

-0.115 c 
-0.284 c 
-1.06 d 
-1.96 d 



-0.6 9 



"Reference pi 
^Reference [122 



e our own KS calculations 
•^Reference [f25| 



c Reference p3 | 
^Reference 



^Reference 26 
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TABLE II. The vibrational frequencies (in cm *) of Na, Mg, Li dimer using the OFM 
and KS self consistent method. 



dimers OFM KS 

Na 167.4 168 

Mg 107.3 108.6 

Li 273.7 311 



TABLE III. The bondlengths between Li-AL of Li n Al,n =1,8 using OFM compared 



with those obtained by KS method [ 20] . All the bondlengths are in a.u 



system 

LiAl 

Li 2 Al 

Li 3 Al 

Li 4 Al 

Li 5 Al 

LiqAI 
Li 7 Al 



Li 8 Al 



OFM 
4.77 

2 x 4.76 

3 x 4.79 

4 x 4.84 

4 x 4.84 
4.95 

6 x 4.79 
4.97 

2 x 4.92 
2 x 4.88 
2 x 4.89 
8 x 4.99 



KS 
5.35 

2 x 5.22 

3 x 4.98 
2 x 4.82 
2 x 4.89 

4 x 4.74 
5.13 

6 x 4.58 
4.70 

2 x 4.85 
2 x 4.74 
2 x 4.81 
8 x 4.82 



% error 

10.8 

8.8 

3.8 

0.4 

1 

2 

3.3 
4.5 
5.7 
1.4 
2.9 
1.6 
3.5 



Symmetry 

n 

^oov 

c 2v 



C 



o h 

Cih 



D 



id 
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Figure Captions 

1. The self-consistent charge densities of Al dimer. Curve a represents the F(N e ) func- 
tional charge density, curve b represents the A = | charge density and curve c denotes 
the charge density obtained using the KS method. 

2. Comparison of self-consistent charge densities by the F(N e ) (curve a) and A = | 
(curve c) functional for Al dimer with the superposed ( th iteration) free Al atom 
charge density (curve b). 

3. The ground state geometries of the Li n Al clusters for n — 1,8. The large sphere 
represents the Li atom and small sphere represents the Al atom. 
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